Cite this Article
Liu Ruo-Yun, Li Zheng-Lin. Effects of rough surface on sound propagation in shallow water. Chinese Physics B, 2019, 28(1): 014302  
Permissions
Effects of rough surface on sound propagation in shallow water
Liu Ruo-Yun1, 2, Li Zheng-Lin1, †
State Key Laboratory of Acoustics, Institute of Acoustics, Chinese Academy of Sciences, Beijing 100190, China
University of Chinese Academy of Sciences, Beijing 100190, China

 

† Corresponding author. E-mail: lzhl@mail.ioa.ac.cn

Project supported by the National Natural Science Foundation of China (Grant Nos. 11434012, 11874061, and 41561144006).

Abstract

Underwater acoustic applications depend critically on the prediction of sound propagation, which can be significantly affected by a rough surface, especially in shallow water. This paper aims to investigate how randomly fluctuating surface influences transmission loss (TL) in shallow water. The one-dimension wind-wave spectrum, Monterey–Miami parabolic equation (MMPE) model, Monte Carlo method, and parallel computing technology are combined to investigate the effects of different sea states on sound propagation. It is shown that TL distribution properties are related to the wind speed, frequency, range, and sound speed profile. In a homogenous waveguide, with wind speed increasing, the TLs are greater and more dispersive. For a negative thermocline waveguide, when the source is above the thermocline and the receiver is below that, the effects of the rough surface are the same and more significant. When the source and receiver are both below the thermocline, the TL distributions are nearly the same for different wind speeds. The mechanism of the different TL distribution properties in the thermocline environment is explained by using ray theory. In conclusion, the statistical characteristics of TL are affected by the relative roughness of the surface, the interaction strength of the sound field with the surface, and the changes of propagating angle due to refraction.

PACS: 43.30.+m;43.30.Hw
Keyword:surface fluctuation;shallow water;transmission loss;statistical characteristics
1. Introduction

Underwater acoustic applications depend critically on the prediction of sound propagation, which can be significantly affected by the wind-caused fluctuating surface, especially in shallow water. Thus, fluctuating surface has important influence on underwater applications, such as underwater acoustic communication,[1,2] inversion of bottom absorption coefficients,[3] and matching-field localization.[4] This implies that it is essential to consider the effect of surface roughness on the performance evaluation of the underwater communication and estimate the influence of the wind-caused sea wave on the experiments before further processing and analyzing the received signals. Therefore, it is quite important to investigate how the randomly fluctuating surface affects sound propagation in shallow water.

Sound propagation is determined by boundary conditions to some extent. Generally, sea surface is modeled as a smooth interface, where the sound is specularly reflected. However, in the case of strong wind, the roughness and variation of the surface should be considered. Diffuse reflection from the rough surface makes sound ray stochastically scattered to disparate directions besides the specular reflection direction, leading to the attenuation and fluctuation of acoustic pressure at the receiver. Other effects include arrival time fluctuation,[5,6] signal coherence varyiation,[7–9] focusing and defocusing,[10] angle spreading,[11] frequency shifting, and bandwidth broadening.[12] In shallow water, since repeated surface–bottom interactions occur, these effects will be more significant.

Over the recent decades, there have been many papers discussing the effects of wind-generated fluctuating surfaces. Early research had reported the backscattering of sound signals for monostatic sonar systems.[13,14] While recent studies focused on forward scattering with the development of bistatic sonar systems. Kuperman found that the attenuation caused by surface scattering was affected by wind speed, sound speed profile, bottom properties, and wind directions.[15] Weston systematically examined the wind–wave–acoustics relation and concluded the fitting formula among scattering attenuation, wind speed, and frequency.[16] Wang reported that the significant difference of transmission loss (TL) in two underwater experiments resulted from the wind-caused rough surface.[3] Using experiments and large numbers of simulations in estuarine shallow water, Zou indicated a strong relation between surface roughness and pressure amplitude of broadband signals at the once-surface-reflected path.[17] Since wind-caused surface fluctuation is a random process, it is effective and reasonable to use statistical methods to study its effects, which will bring the challenge of large computing requirement in the simulation. Therefore, few researchers use these methods. While we have employed parallel computing to obtain the simulation results to be analyzed with these methods.

This paper uses a one-dimensional wind-wave spectrum for fluctuating surface generation and Monterey–Miami parabolic equation (MMPE) model to simulate the sound propagation under the rough surface. A parallel algorithm is proposed to speed up the simulation process and the Monte Carlo method is used to model the stochastic sound field. The purpose of this paper is to discuss the statistic characters of TL distribution caused by rough surface scattering under different wind-wave conditions in shallow water, and to introduce a rapid approach to evaluate the variation of TL for the fluctuating surface. This study is based on numerical simulations of an isothermal environment and a thermocline environment. As an innovation, we propose a rapid method to calculate the statistical properties of TL under randomly fluctuating surface.

The remainder of the article is organized as follows. Section 2 introduces the combined numerical model based on a one-dimensional wind-wave spectrum, MMPE model, Monte Carlo method, and parallel computing technology. Section 3 presents the simulation results of fluctuation of sound fields caused by a rough surface for different wind speeds. And the influencing factors of the different TL distribution properties are discussed in Section 3. The conclusions are summarized in Section 4.

2. Method
2.1. Generation of rough surface

In this paper, a one-dimensional wind-wave spectrum is used to generate stochastically fluctuating surfaces.[18] In two-dimensional coordinates (r,z), the surface interface is defined as follow:

where η(r,t) is the surface displacement, a random function of range and time, increasing downward. Using the following method to generate the displacement function η(r) at a certain moment, where η(r) is defined as a zero-mean Gaussian random function with
Here k is the spatial wavenumber. W1(k) is the roughness spectrum, and its approximate formation is
where β is the weighted index of the spectrum and Lcorr is the correlation length
with σ being the root mean square (RMS) of the wave height of the rough surface. Assuming the direction of the wind is along the coordinate r, an empirical formula of the relation between σ and wind speed v at 20 m above the sea surface[19] is
where g is the acceleration of gravity.

In order to compute stochastic realizations of roughness, superimpose a component Aeiθ upon the wave spectrum W1(k), consisting of a random wavenumber amplitude and phase θ = 2 πR1, where R1 and R2 are uniformly distributed random variables within the interval (0,1). And this complex spectrum is conversed into r-domain with discrete fast Fourier transform (FFT) to yield the surface displacement function η(r),

And then η(r) is normalized by its actual RMS value, i.e., η(r)/⟨η2(r)⟩1/2, and multiplied by the desired RMS value σ. The result of this calculation is a single stochastic realization of a rough surface η(r). For example, a randomly fluctuating surface with wind speed of 15 m/s and weighted index of the spectrum β = 3.5 is generated according to this approach, as displayed in Fig. 1.

Fig. 1. A randomly realization of the fluctuating surface for wind speed of 15 m/s.
2.2. MMPE model with rough surfaces

The MMPE acoustic propagation model was proposed by Kevin B. Smith[20] based on UMPE model[18] using the split-step Fourier (SSF) algorithm[21] to solve the parabolic equation (PE). The MMPE model is only capable of computing fields scattered from a flat surface. So, a field transformation approach that allows for the calculation of rough surface scattering is used in this paper, which was introduced by Tappert and Nghiem-Phu.[22]

In the MMPE Model, the waveguide is considered as a constant-density medium in cylindrical coordinates (r,z,φ) and a point source of unit strength is set at coordinates (r = 0, z = zs, φ = 0). The Helmholtz equation is

where p is the sound pressure, k0 = ω/c0 is the reference wavenumber (ω is the angular frequency), and n(r,z,φ) = c0/c(r,z,φ) is the index of refraction. c(r,z,φ) is the sound speed and c0 is the reference sound speed. Assume the solution of Eq. (7) is in the form of an outward cylindrical wave
where u(r,z) is the pressure function in two dimensions and the term accounts for azimuthal spreading. Substituting Eq. (8) into the two-dimensions Helmholtz equation, which is obtained by ignoring the source term and the azimuthal coupling term in Eq. (7), the parabolic equation can be derived as
which represents the forward propagating sound energy in the waveguide. ψ(r,z) is called the PE field function, or envelope function, and is defined as
The approximation of the operator Qop employed in the MMPE model was introduced by Thomson and Chapman,[23]
It is a higher order approximation than the standard PE, called the wide-angle approximation (WAPE). Equation (9) can be solved by giving a source-field distribution over depth at the initial range and hence utilizing a range-marching numerical technique, split-step Fourier (SSF) algorithm. The SSF solution of the parabolic equation is
where
Uop (r,z) is in z-domain and is in kz-domain (wavenumber domain). The convention of the two domains is employed when using FFT. The starting field for the wide-angle point source is given by
where R0 = 1 m, and zs is the source depth.

To satisfy the boundary condition of smooth pressure release surface, ψ(z = 0) = 0, assume an identical image ocean overlays the real ocean for the negative values of depth. The PE field function is exactly equal but of opposite sign in the image ocean, i.e.,

While in the modified approach for the rough surface forward scattering,[22] the field function is an odd function about the rough surface interface z = η(r). So, for the modified PE field function, and
where ψ(r,z) is computed with Eq. (12), and the generation approach of η(r) is introduced in the previous section. And the expression of operator Uop is rewritten as
According to Eqs. (8) and (10), the final pressure expression is

2.3. Monte Carlo method

The Monte Carlo method, also called the statistical simulation method, is to use statistic experiment of random variables to simulate the solution of mathematic, physical, and engineering problems.

For the randomness of the sea surface, from which the sound wave is scattered, the TL at the receiver is randomly fluctuating as well. The Monte Carlo method is used to obtain the distribution of TL. According to Subsection 2.1, rough surfaces are generated in 100 different stochastic realizations based on two random variables R1 and R2 to model the varying surface displacement functions η(r). Then the sound field under each one of the 100 surfaces is calculated separately according to Subsection 2.2, which is used to simulate the different pressures at the receiver for different stochastic realizations. Next, the TL is evaluated, and the distribution is calculated.

Generally, the TL of a single frequency signal is vibrating with the range because of the interference between normal modes. To get smoother TL curve, we average the sound energy of N discrete frequencies within 1/3 octave band width around the center frequency,

where pi (r,z) is the pressure of each frequency, and N = 10 in this paper. Thus, each Monte Carlo experiment consists of 1000 field calculations, and a long processing time is expected.

2.4. A rapid parallel algorithm

To reduce the processing time, a rapid algorithm is proposed based on parallel computing, in which multiple processes are executed simultaneously to complete a large number of calculations quickly. The parallel algorithm is based on a message passing interface (MPI), a message passing library specification. The library routines are called from Fortran programs to accomplish communication between processes.[24]

In this algorithm, a subtask is defined as sound field prediction of a single frequency under a rough surface. Executed on m cores of a multi-core computer, it is most time saving to allocate all subtasks equally to the m processes. The field calculating procedure of each subtask is independent of the others. Thus, the problem could be easily parallelized and communication between processes only occurs at the start and the end of the program to deliver the parameters and gather the results.

The parallel MMPE algorithm is designed in the master–slave mode, as shown in Fig. 2. Both master and slave processes perform the sound field calculation. The master process is responsible for operations, including extra managing processes, assigning tasks, carrying out file operations, broadcasting the parameters to the others, and gathering the results.

Fig. 2. The parallel MMPE program on a computer with multiple cores.

The validity of the parallel MMPE program is illustrated through a comparison with the RAM model. In the waveguide of Fig. 5, a 1000 Hz source is set at 7 m below the sea level and a receiver is set at 10 m below the sea level. The TL is displayed in Fig. 3 for an environment with a flat surface, which shows that the results agree well.

Fig. 3. Comparison of the TL results from MMPE and RAM.
Fig. 4. (a) Speedup and (b) parallel efficiency of the MMPE parallel program.
Fig. 5. Shallow water environment in the South China Sea.

The performance of the parallel program is measured by both the speedup ratio Sm and the parallel efficiency Em, which are defined as follows:

where T1 and Tm are the computing time for the serial program and for the parallel program executed on m processors, respectively. Figure 4 shows the ideal and actual speedup and parallel efficiency on a four-core computer with m equal to 1, 2, 3, and 4, separately. The processing time of the parallel part reduces in the multi-processor case. In ideal situation, Sm is equal to m and Em is always 100%. While the processing time of the communication part and serial part, such as file operations, is not less in the multiple processors case than that in the one processor case. For the existence of these parts in the algorithm, the actual speedup and parallel efficiency are always less than the ideal ones. In Fig. 4, as the number of processors rises, the actual speedup increases nearly linearly, close to the ideal case, and the parallel efficiency declines slightly. The best speedup available on this computer is up to 3.47 in the 4-processor case, meanwhile the parallel efficiency is 87%. It is shown that the performance of the parallel program is satisfactory and the processing time is significantly less than that of the serial program.

3. Fluctuation of sound fields caused by rough surface for different wind speeds
3.1. Simulation results in an isothermal environment

A typical range-independent homogenous waveguide in the South China Sea is used,[3] as displayed in Fig. 5. The sound speed of the water layer is 1537 m/s, the depth is 88 m, and the attenuation coefficient of water α1 is calculated using an empirical expression dependent on frequency (f in kHz)

The bottom consists of a sediment layer with thickness of 8 m, constant compressional sound speed of 1606 m/s, density of 1.65 g/cm3, and compressional attenuation of 0.517(f/1000)1.07 dB/λ, and a homogeneous half-space basement with sound speed of 1704 m/s, density of 1.9 g/cm3, and the same attenuation of the sediment. The source depth is 7 m, and the receiver depth is 10 m. Using Eq. (21), TLs are calculated at 200 Hz and 1000 Hz respectively. Two source–receiver ranges of 20 km and 50 km are selected.

To investigate the regulation of TL distribution under the rough surface, simulation is conducted for four different wind speeds (v = 7.5 m/s, 10 m/s, 15 m/s, and 20 m/s) respectively, corresponding to four, five, seven, and eight wind scales. A set of 100 randomly fluctuating surfaces of each wind speed are separately generated by the wind-wave spectrum. The sound field under each set of surfaces is calculated and the result is obtained by averaging the energy of 10 frequencies within 1/3 octave bandwidth, according to Eq. (20). Then the TLs of each set are divided into several groups, of which the interval length is 1 dB. Assume TL is the interval center of a certain group, and then the probability q(TL) is calculated with the following equation:

where N1(TL) is the number of samples in this group and N0 is the number of total samples. Figure 6 shows the probability of TL for 200 Hz when the wind speed is 20 m/s. TL probability curves under different wind speeds are shown in Fig. 7.

Fig. 6. Probability of TL at 200 Hz with wind speed of 20 m/s.
Fig. 7. Probability of TL at different frequencies and different ranges for the receiver at 10 m, where panels (a) and (b) are for frequency 200 Hz at the ranges of 20 km and 50 km, respectively; panels (c) and (d) are for frequency 1000 Hz at the ranges of 20 km and 50 km, respectively.

In Fig. 6, the TL of each range randomly distributes in an interval of 10 dB to 15 dB. The interval length increases from 0 to 20 km, for the randomness increase with the number of fluctuating surface reflections. The maximal TL probability changes from 0.18 to 0.3 in 60 km, The TL near the edge has less probability, less than 0.05.

In Fig. 7, the TL of each case distributes in an interval around a maximum point corresponding to the peak probability. We define the interval length as the variation range of TLs. We can use the interval length and peak probability to evaluate the dispersion degree and use the maximum point to measure the magnitude of TL. Obviously, the TL is increasing and more dispersive as the wind speed rises. For instance, in the case of 200 Hz at 50 km, displayed in Fig. 7(b), the peaks of TL probability for four wind speeds are at 78.5 dB, 79.5 dB, 82.5 dB, and 87.5 dB, respectively, which represents the rise of TL. Because higher wind speed leads to larger RMS wave height according to Eq. (5). And a rougher surface increases the scattering attenuation. Besides, the distribution of the TL for the four wind speeds varies within an interval of 2 dB, 5 dB, 12 dB, and 14 dB, respectively. And the peak probabilities are respectively 0.74, 0.58, 0.25, and 0.26. The interval increases and the peak probability tends to decline, which indicates that the TL distribution is more dispersive in the case of stronger wind. As the surface varies more dramatically for different stochastic realizations, the changing of surface scattering paths is more complex, making the sound field more fluctuating.

Compared with the case of 200 Hz at 20 km in Fig. 7(a), the TL increase with wind speed is more significant for the TL at 50 km in Fig. 7(b). It is the same for the 1000 Hz case in Figs. 7(c) and 7(d). For example, the maximum points of 7.5 m/s and 20 m/s differ by 6 dB in Fig. 7(a), while 9 dB in Fig. 7(b). It is due to sound wave going through more surface reflections before arriving at 50 km. Thus, the scattering attenuation is more significant.

From the comparison between Figs. 7(a) and 7(c) (or Figs. 7(b) and 7(d)), we can find that as the wind speed rises, the increase of scattering attenuation is more significant in the 1000 Hz case than that in 200 Hz. For instance, when the wind speed rises from 7.5 m/s to 20 m/s, the maximum points of TL increase 9 dB for 200 Hz and 11 dB for 1000 Hz at 50 km. This implies that the sound wave with higher frequency decays faster under the same rough surfaces. In acoustics, the scale of length is relative to the sound wavelength. For the same absolute value of surface displacement, the relative value is larger in the 1000 Hz case than that in the 200 Hz case. Thus, the effect of rough surface is more significant for higher frequency.

Besides the TL caused by rough surface TLs, the TL caused by volume attenuation TLw also increases with frequency. The TL of the receiver depth, 10 m, is use to evaluate TLs

where TLrough(f,v) is the TL for the frequency of f under the wind speed of v, and TLflat (f) is under no-wind conditions. And TLw is calculated with the following equation
where α1 is the attenuation coefficient of water in Eq. (23). At the range of 50 km, TLs is 8.5 dB and TLw is 0.40 dB for 200 Hz. And TLs is 11 dB and TLw is 3.47 dB for 1000 Hz. TLs increases by 2.5 dB and TLw increases by 3.07 dB. The comparison implies that for 200 Hz and 1000 Hz, rough surface is the main attenuation factor, while with the frequency rising, TLs and TLw increase by the same level.

The simulation results in this isothermal environment indicate that the statistical characteristics of TL is affected by the wind speed, range, and frequency.

3.2. Simulation results in a thermocline environment

There is a typical waveguide with negative thermocline layer from the Yellow Sea experiment in 1996,[25] of which the environment parameters are shown in Fig. 8. The sound speed above the thermocline is 1537 m/s, and the underneath is 1480 m/s. Within the thermocline layer from 13 m to 28 m, a linear sound-speed profile with a constant negative gradient of −3.8 s−1 is used. The attenuation coefficient α1 in water is shown in Eq. (23). The central frequency is 1000 Hz. Two source depths are set at 7 m and 50 m respectively to represent above and below the thermocline. And the receiver is placed at 70 m, below the thermocline, and the source–receiver range is 20 km. The TL probability curves for four wind speeds (v = 7.5 m/s, 10 m/s, 15 m/s, and 20 m/s) are shown in Fig. 9.

Fig. 8. Shallow water environment with a thermocline.
Fig. 9. Probability of TL in a thermocline environment at 1000 Hz, where panel (a) is for source above the thermocline (Zs = 7 m, ZR = 70 m), and panel (b) is for source below the thermocline (Zs = 50 m, ZR = 70 m).

Figure 9(a) shows the statistical characteristics for the case of a source located above the thermocline and a receiver located below the thermocline. The TLs are larger and more dispersive as the wind speed rises, the same as the homogeneous waveguide. Compared with Fig. 7(c), the effects of rough surface in this case are more significant than the homogeneous waveguide. For example, the increment of maximum points is 6 dB from 7.5 m/s to 20 m/s in Fig. 7(c), while it is 21 dB in Fig. 9(a), which means that in this case the increase of TL is greater with wind speed rising. And under 20 m/s, the peak of probability is 0.23 and the variation range of TLs is 13 dB in Fig. 7(c), and they are 0.17 and 17 dB in Fig. 9(a), which demonstrates that the distribution of TL is more dispersive in this case. When the source is above the thermocline, sound rays propagate in the whole depth, reflected by the surface and bottom, which are called surface-reflected bottom-reflected (SRBR) rays. Because of the bottom attenuation, the rays with greater angles decay quickly after several bottom reflections. For the environment shown in Fig. 8, the rays with grazing angles less than the critical angle 21.16° can propagate long distances. In Fig. 10, some typical rays with small grazing angles are plotted for wind speeds of 0 m/s, 7.5 m/s, and 20 m/s, respectively. In Fig. 10(a), these rays can reach 20 km in the environment with a flat surface. In Figs. 10(b) and 10(c), the reflection angles of these rays are changed after the reflection of the rough surface. Those rays with increased reflection angles lead to large bottom loss. Thus, they cannot propagate to 20 km, and the TL of 20 km increases. For the wind speed of 20 m/s, more rays vanish before arriving at 20 km. Besides, lots of rays are bent downward in the thermocline, and then they incident into the bottom with greater grazing angles. The thermocline enlarges the effects of rough surface. Therefore, the increase of TL in the thermocline environment is greater than that in the homogenous waveguide, and TL is more dispersive.

Fig. 10. Rays for source above the thermocline (Zs = 7 m) with wind speeds of (a) 0 m/s, (b) 7.5 m/s, and (c) 20 m/s.

In Fig. 9(b), with the source and receiver located beneath the negative thermocline, the TL is quite concentrated within an interval of 1 dB for all wind speeds. Figure 11 shows the eigen rays to the receiver for wind speeds of 7.5 m/s and 20 m/s, respectively. Most of the sound rays propagating upward are refracted and bent downward to the bottom in the negative thermocline layer without interaction with the surface, called refracted bottom-reflected (RBR) rays, such as the green, red, and blue rays in Fig. 11. The rest is SRBR rays with grazing angle greater than 15.65°, like the black ray in Fig. 11(a). They are affected by the rough surface, attenuating quickly and having less energy. The acoustic rays arriving at 70 m mainly consist of RBR rays and few SRBR rays. Therefore, the total sound energy is rarely affected by fluctuating surfaces, and then the TL distributions for different wind speeds are almost the same.

Fig. 11. Eigen rays for source below the thermocline (Zs = 50 m, ZR = 70 m) with the wind speeds of (a) 7.5 m/s and (b) 20 m/s.
4. Conclusion

Under a randomly fluctuating surface, a sound field in shallow water attenuates and randomly distributes in an interval. The simulations indicate that the statistical characteristics of TL are affected by the wind speed, range, frequency, and sound speed profile. Source and receiver depths relative to the thermocline are influential as well.

The TL is increasing and more dispersive as the wind speed raises due to the roughness of the surface. The increase is more significant for the higher frequency case due to the relative roughness of the surface. The interaction strength with the surface accounts for the following two phenomena. The TL increase with wind speed is more significant for the receiver farther away from the source. For the case of the source below the thermocline, the TL for the receiver below the thermocline is barely influenced by rough surfaces. While the TL for the source above the thermocline is more affected than that in an isothermal environment, because the grazing angle increases in the thermocline. In conclusion, the influencing factors are the relative roughness of the surface, the interaction strength of the sound field with the surface, and the changes of propagating angle due to refraction.

Besides, a reliable and rapid approach is introduced to simulate the distribution of TL under randomly rough surface environment. It can be used to evaluate the average value and the variation of TL under different wind-wave conditions. The regulations and the approach discussed in this paper is expected to be helpful for other underwater acoustic studies, such as sound field fluctuation with the presence of internal waves.

Reference
[1] Luo J Badiey M 2005 J. Acoust. Soc. Am. 118 2038
[2] Al-Kurd A 1996 J. Acoust. Soc. Am. 100 2703
[3] Wang X H Peng Z H Li Z L 2007 Tech. Acoust. 26 551 in Chinese
[4] Li Z L 2001 The Effect of Internals Waves, Surface Fluctuation and Bottom Roughness on Matched Field Source Localization in Shallow Water Ph. D. dissertation Beijing Graduate University of Chinese Academy of Sciences in Chinese
[5] Park C Seong W Gerstoft P Hodgkiss W S 2011 J. Acoust. Soc. Am. 129 98
[6] Karjadi E A Badiey M Kirby J T Bayindir C 2012 IEEE J. Ocean. Eng. 37 112
[7] Badiey M Mu Y K Simmen J A Forsythe S E 2000 IEEE J. Ocean. Eng. 25 492
[8] Yang T C 2006 J. Acoust. Soc. Am. 120 2595
[9] Dahl P H 1996 J. Acoust. Soc. Am. 100 748
[10] Tindle C T Deane G B 2005 J. Acoust. Soc. Am. 117 2783
[11] Dahl P H 2001 IEEE J. Ocean. Eng. 26 141
[12] Siderius M Porter M B 2008 J. Acoust. Soc. Am. 124 137
[13] Eckart C 1953 J. Acoust. Soc. Am. 25 566
[14] Urick R Hoover R 1956 J. Acoust. Soc. Am. 28 1038
[15] Kuperman W A Ingenito F 1977 J. Acoust. Soc. Am. 61 1178
[16] Weston D E Ching P A 1989 J. Acoust. Soc. Am. 86 1530
[17] Zou Z G Badiey M 2018 IEEE J. Ocean. Eng. 43 1187
[18] Smith K B Tappert F D 1993 UMPE: University Miami Parabolic Equation Model Version 1.0
[19] Neumann G Pierson W J 1966 Principle of Physical Oceanography Englewood Cliffs Prentice Hall 351
[20] Smith K B 2001 J. Comput. Acoust. 9 243
[21] Hardin R H Tappert F D 1973 SIAM Rev. 15 423
[22] Tappert F D Lan N 1985 J. Acoust. Soc. Am. 77 S101
[23] Thomson D J Chapman N R 1983 J. Acoust. Soc. Am. 74 1848
[24] Lusk E Gropp W 1995 Adv. Parallel Computing 10 265
[25] Li F H Zhang R H 2000 Acta Acust. 25 297 in Chinese